clear;clc;
A=[0.98 0.1 0;0.02 0.7 0.05;0 0.2 0.95];
T=1000000;
S=zeros(1,T);
S(1)=2;
[b,v]=eig(A);
for t=1:T-1
   S(t+1)=myfun(S(t),A);
end
figure;
hist(S);
a1=length(find(S==1));
a2=length(find(S==2));
a3=length(find(S==3));
A1=1;A2=a2/a1;A3=a3/a1;
B1=1;B2=b(2,2)/b(1,2);B3=b(3,2)/b(1,2);
figure;plot([1,2,3],[A1 A2 A3],'-*',[1 2 3],[B1 B2 B3],'-+');